Testing Out Hidden Markov Models

Establishing a mixture distribution made up of two gaussians

(well separated means)

X <- rnorm(1000, 1, 0.5)
Y <- rnorm(1000, 10, 2)
Z <- rbind(X,Y)
hist(Z, breaks=30, xlab="Value")

Establish two-state time course data

set.seed(seed=2)
X1 <- rnorm(10000, 1, 0.3)
X2 <- rnorm(10000, 10, 0.6)
trans_probs <- matrix(data=c(0.95, 0.05, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
data <- data.frame(X1, X2)

n = 10000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj <- sapply(1:n, function(x) data[x, t_states[x]])

plot(traj[1:1000], pch=16, ylab="Value", xlab="Time")

Model the dataset using a two state Hidden Markov Model

hmm <- depmix(response = traj ~ 1, data=data.frame(traj), nstates=2, instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: -26440.71 
## iteration 5 logLik: -16554.2 
## converged at iteration 10 with logLik: -8730.31

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.955 0.045
## fromS2 0.009 0.991
## 
## Response parameters 
## Resp 1 : gaussian 
##     Re1.(Intercept) Re1.sd
## St1           1.013  0.300
## St2           9.998  0.606

esttrans <- posterior(f)

Visualize Results

plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', ylab=c("Value"), xlab=c("Time"))

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', ylab=c("Value"), xlab=c("Time"))

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", ylab=c("Value"), xlab=c("Time"))

New mixture with closer means

X <- rnorm(10000, 1, 0.5)
Y <- rnorm(10000, 2, 0.6)
Z <- rbind(X,Y)
hist(Z, breaks=30, xlab="Value")

Establish two-state time course data

##Establish Two-State Time Course Data
set.seed(seed=3)
X1 <- rnorm(10000, 1, 0.5)
X2 <- rnorm(10000, 2, 0.5)
trans_probs <- matrix(data=c(0.95, 0.05, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
data <- data.frame(X1, X2)

n = 10000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj <- sapply(1:n, function(x) data[x, t_states[x]])

plot(traj[1:1000], pch=16, xlab="Time", ylab="Value")

Model the dataset using a two state Hidden Markov Model

hmm <- depmix(response = traj ~ 1, data=data.frame(traj), nstates=2, instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: -9332.465 
## iteration 5 logLik: -9323.201 
## iteration 10 logLik: -9258.26 
## iteration 15 logLik: -9083.684 
## iteration 20 logLik: -8721.498 
## iteration 25 logLik: -8105.488 
## iteration 30 logLik: -7962.52 
## converged at iteration 34 with logLik: -7962.476

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.991 0.009
## fromS2 0.056 0.944
## 
## Response parameters 
## Resp 1 : gaussian 
##     Re1.(Intercept) Re1.sd
## St1           1.989  0.508
## St2           0.998  0.509

esttrans <- posterior(f)

Visualize Results

plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', xlab="Time", ylab="Value")

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Time", ylab="State")

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Time", ylab="State")

Construct a hypothetical dataset using microscopy data

##Use 1D measures so that they are gaussian
x <- Sox2_data_unsync_20sec
L_fitX <- fitdistr(x[x$Cell_Line == "121316.2E1_ESC", "Corrected X_Distance (um)"], densfun = "normal")
L_fitY <- fitdistr(x[x$Cell_Line == "121316.2E1_ESC", "Corrected Y_Distance (um)"], densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Corrected X_Distance (um)"], freq = FALSE, ylim=c(0,5), xlab="X Distance (um)", xlim=c(-0.5,0.5), main="Looped Model X Distance Distribution")
L_model <- dnorm(step, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)

hist(x[x$Cell_Line == "121316.2E1_ESC", "Corrected Y_Distance (um)"], freq = FALSE, ylim=c(0,5), xlab="Y Distance (um)", xlim=c(-0.5,0.5), main="Looped Model Y Distance Distribution")
L_model <- dnorm(step, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)

U_fitX <- fitdistr(x[x$Cell_Line == "013116.2G8_ESC", "Corrected X_Distance (um)"], densfun = "normal")
U_fitY <- fitdistr(x[x$Cell_Line == "013116.2G8_ESC", "Corrected Y_Distance (um)"], densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "013116.2G8_ESC", "Corrected X_Distance (um)"], freq = FALSE, ylim=c(0,2.5), xlab="X Distance (um)", xlim=c(-0.75,0.75), main="Unlooped Model X Distance Distribution")
U_model <- dnorm(step, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")

hist(x[x$Cell_Line == "013116.2G8_ESC", "Corrected Y_Distance (um)"], freq = FALSE, ylim=c(0,2.5), xlab="Y Distance (um)", xlim=c(-0.75,0.75), main="Unlooped Model Y Distance Distribution")
U_model <- dnorm(step, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")

##Establish as hypothetical two-state time course based on the data

set.seed(seed=4)
L_X <- rnorm(10000, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
U_X <- rnorm(10000, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])

L_Y <- rnorm(10000, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
U_Y <- rnorm(10000, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])

trans_probs <- matrix(data=c(0.95, 0.01, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
Xdata <- data.frame(L_X, U_X)
Ydata <- data.frame(L_Y, U_Y)

n = 10000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])

plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Distance (um)")

plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Distance (um)")

Model the dataset using a two state Hidden Markov Model

hmm <- depmix(response = list(traj_X ~ 1, traj_Y ~ 1), ntimes=10000, nstates=2, family=list(gaussian(), gaussian()), instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: 4044.465 
## iteration 5 logLik: 4211.006 
## iteration 10 logLik: 4876.121 
## iteration 15 logLik: 5705.837 
## iteration 20 logLik: 6055.305 
## converged at iteration 25 with logLik: 6055.491

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.989 0.011
## fromS2 0.012 0.988
## 
## Response parameters 
## Resp 1 : gaussian 
## Resp 2 : gaussian 
##     Re1.(Intercept) Re1.sd Re2.(Intercept) Re2.sd
## St1           0.015  0.121           0.016  0.126
## St2           0.034  0.258          -0.012  0.246

esttrans <- posterior(f)

Visualize results

plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")

Function for generating autocorrelated time series

autoFunc <- function(x, rho, dist){
  
  y <- rho * x + sqrt(1 - rho^2) * sample(dist, size=1)
  return(y)
}

set.seed(5)
z1 <- rnorm(n=1000, mean=1.5, sd=2)
X <- array(dim=length(z1))
cor <- 0.8
X[1] <- z1[1]
for (i in 2:length(z1)){
  X[i] <- autoFunc(X[i - 1], cor, z1)
}

z2 <- rnorm(n=1000, mean=3, sd=2)
Y <- array(dim=length(z2))
cor <- 0.8
Y[1] <- z2[1]
for (i in 2:length(z2)){
  Y[i] <- autoFunc(Y[i - 1], cor, z2)
}

plot(X, pch=16, col="red", ylim=c(0, 15), ylab="Value", xlab="Time", main="Two autocorrelated time course datasets from which to sample")
legend(x="topright", pch=16, legend=c("State1", "State2"), col=c("red", "purple"))
points(Y, pch=16, col="purple")

Setup mixture dataset

trans_probs <- matrix(data=c(0.99, 0.01, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
data <- data.frame(X, Y)

n = 10000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj <- sapply(1:n, function(x) data[x, t_states[x]])

plot(traj[1:1000], pch=16, xlab="Time", ylab="Value", main="Autocorrelated time course dataset derived from two independent states")

Model the dataset using a two state Hidden Markov Model

hmm <- depmix(response = traj~1, ntimes=10000, nstates=2, family=gaussian(), instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: -2516.253 
## iteration 5 logLik: -2516.207 
## iteration 10 logLik: -2516.108 
## iteration 15 logLik: -2515.786 
## iteration 20 logLik: -2514.695 
## iteration 25 logLik: -2510.865 
## iteration 30 logLik: -2497.06 
## iteration 35 logLik: -2440.09 
## iteration 40 logLik: -2347.081 
## iteration 45 logLik: -2279.898 
## iteration 50 logLik: -2235.285 
## iteration 55 logLik: -2206.705 
## iteration 60 logLik: -2188.935 
## iteration 65 logLik: -2178.302 
## iteration 70 logLik: -2172.172 
## iteration 75 logLik: -2168.757 
## iteration 80 logLik: -2166.91 
## iteration 85 logLik: -2165.935 
## iteration 90 logLik: -2165.43 
## iteration 95 logLik: -2165.172 
## iteration 100 logLik: -2165.041 
## iteration 105 logLik: -2164.975 
## iteration 110 logLik: -2164.943 
## iteration 115 logLik: -2164.926 
## iteration 120 logLik: -2164.918 
## iteration 125 logLik: -2164.914 
## iteration 130 logLik: -2164.912 
## iteration 135 logLik: -2164.911 
## iteration 140 logLik: -2164.91 
## iteration 145 logLik: -2164.91 
## converged at iteration 149 with logLik: -2164.91

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.974 0.026
## fromS2 0.037 0.963
## 
## Response parameters 
## Resp 1 : gaussian 
##     Re1.(Intercept) Re1.sd
## St1           4.365  1.873
## St2           9.031  2.023

esttrans <- posterior(f)

Visualize results

plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', xlab="Time", ylab="Value")

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Time", ylab="State")

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Time", ylab="State")

Conclusion: HMM approach fails when data has autocorrelation.

Using Displacements of Data

Generate Gaussian Displacments

x <- Sox2_data_unsync_20sec
images <- unique(x$C1_Image)

for (i in 1:length(images)) {
  start = 1;
  end = 100;
  cur_image = images[i];
  cur_data = x[x$C1_Image == cur_image,];
  for (j in start:end) {
    if (nrow(cur_data[cur_data$`C1_T Value (frame)` == (j - 1),]) != 0 & nrow(cur_data[cur_data$`C1_T Value (frame)` == j,]) != 0){
      x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_X.Displacement..um."] <- (x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == (j - 1), "Relative_C1_Corrected_X Value (um)"] - x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_C1_Corrected_X Value (um)"]);
      x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_Y.Displacement..um."] <- (x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == (j - 1), "Relative_C1_Corrected_Y Value (um)"] - x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_C1_Corrected_Y Value (um)"]);
    }
    else{
      x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_X.Displacement..um."] <- NA;
      x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_Y.Displacement..um."] <- NA;
    }
  }
}

hist(x[x$Cell_Line == "101515.C3_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

hist(x[x$Cell_Line == "101515.C3_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))

Construct a hypothetical dataset using microscopy data

##Use 1D measures so that they are gaussian
L_fitX <- fitdistr(na.omit(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."]), densfun = "normal" )
L_fitY <- fitdistr(na.omit(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."]), densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."], freq = FALSE, ylim=c(0,5), xlab="X Displacement (um)", xlim=c(-0.5,0.5), main="Looped Model X Displacement Distribution")
L_model <- dnorm(step, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)

hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."], freq = FALSE, ylim=c(0,5), xlab="Y Displacement (um)", xlim=c(-0.5,0.5), main="Looped Model Y Displacement Distribution")
L_model <- dnorm(step, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)

U_fitX <- fitdistr(na.omit(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."]), densfun = "normal")
U_fitY <- fitdistr(na.omit(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."]), densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."], freq = FALSE, ylim=c(0,2.5), xlab="X Displacement (um)", xlim=c(-0.75,0.75), main="Unlooped Model X Displacement Distribution")
U_model <- dnorm(step, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")

hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."], freq = FALSE, ylim=c(0,2.5), xlab="Y Displacement (um)", xlim=c(-0.75,0.75), main="Unlooped Model Y Displacement Distribution")
U_model <- dnorm(step, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")

###Establish as hypothetical two-state time course based on the data

#set.seed(seed=4)
L_X <- rnorm(10000, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
U_X <- rnorm(10000, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])

L_Y <- rnorm(10000, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
U_Y <- rnorm(10000, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])

trans_probs <- matrix(data=c(0.98, 0.02, 0.02, 0.98), byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
Xdata <- data.frame(L_X, U_X)
Ydata <- data.frame(L_Y, U_Y)

n = 10000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])

plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Displacement (um)")

plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Displacement (um)")

Model the dataset using a two state Hidden Markov Model

hmm <- depmix(response = list(traj_X ~ 1, traj_Y ~ 1), ntimes=10000, nstates=2, family=list(gaussian(), gaussian()), instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: 13646.79 
## iteration 5 logLik: 13646.92 
## iteration 10 logLik: 13649.53 
## iteration 15 logLik: 13683.9 
## iteration 20 logLik: 13763.51 
## iteration 25 logLik: 13800.88 
## iteration 30 logLik: 13831.8 
## iteration 35 logLik: 13870.12 
## iteration 40 logLik: 13919.41 
## iteration 45 logLik: 13983.4 
## iteration 50 logLik: 14063.69 
## iteration 55 logLik: 14140.3 
## iteration 60 logLik: 14166.24 
## iteration 65 logLik: 14168.01 
## iteration 70 logLik: 14168.06 
## converged at iteration 73 with logLik: 14168.06

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.977 0.023
## fromS2 0.021 0.979
## 
## Response parameters 
## Resp 1 : gaussian 
## Resp 2 : gaussian 
##     Re1.(Intercept) Re1.sd Re2.(Intercept) Re2.sd
## St1           0.000  0.092          -0.005  0.094
## St2          -0.004  0.138          -0.002  0.150

esttrans <- posterior(f)

Visualize results

plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")

Construct a hypothetical dataset using microscopy data with observed autocorrelation

#set.seed(seed=4)
L_X <- x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."]
U_X <- x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."]

L_Y <- x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."]
U_Y <- x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."]

trans_probs <- matrix(data=c(0.98, 0.02, 0.02, 0.98), byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;

## Load data into DataFrame
Xdata <- data.frame(L_X[1:5000], U_X[1:5000])
Ydata <- data.frame(L_Y[1:5000], U_Y[1:5000])

n = 5000
t_states <- s0;
for (i in 1:(n-1)){
  s <- sample(states, size=1, prob=trans_probs[s, ])
  t_states <- c(t_states, s)
}

traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])

plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Displacement (um)")

plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Displacement (um)")

hmm <- depmix(response = list(traj_X ~ 1, traj_Y ~ 1), ntimes=5000, nstates=2, family=list(gaussian(), gaussian()), instart=c(1,0))
f <- fit(hmm)

## iteration 0 logLik: 6696.141 
## iteration 5 logLik: 7041.577 
## iteration 10 logLik: 7166.998 
## iteration 15 logLik: 7228.651 
## iteration 20 logLik: 7250.898 
## iteration 25 logLik: 7257.632 
## iteration 30 logLik: 7259.593 
## iteration 35 logLik: 7260.171 
## iteration 40 logLik: 7260.346 
## iteration 45 logLik: 7260.4 
## iteration 50 logLik: 7260.418 
## iteration 55 logLik: 7260.423 
## iteration 60 logLik: 7260.425 
## converged at iteration 65 with logLik: 7260.425

summary(f)

## Initial state probabilties model 
## pr1 pr2 
##   1   0 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.954 0.046
## fromS2 0.121 0.879
## 
## Response parameters 
## Resp 1 : gaussian 
## Resp 2 : gaussian 
##     Re1.(Intercept) Re1.sd Re2.(Intercept) Re2.sd
## St1               0  0.089          -0.001  0.090
## St2               0  0.175           0.000  0.178

esttrans <- posterior(f)

Visualize results

plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")

plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")

plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")